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ABSTRACT 

We run numerical simulations of the disruption of satellite galaxies in a Galactic 
potential to build up the entire stellar halo, in order to investigate what the next 
generation of astrometric satellites will reveal by observing the halo of the Milky Way. 
We generate artificial DIVA, FAME and GAIA halo catalogues, in which we look 
for the signatures left by the accreted satellites. We develop a method based on the 
standard Friends-of- Friends algorithm applied to the space of integrals of motion. We 
find this simple method can recover about 50% of the different accretion events, when 
the observational uncertainties expected for GAIA are taken into account, even when 
the exact form of the Galactic potential is unknown. The recovery rate for DIVA 
and FAME is much smaller, but these missions, like GAIA, should be able to test 
the hierarchical formation paradigm on our Galaxy by measuring the amount of halo 
substructure in the form of nearby kinematically cold streams with for example, a 
two-point correlation function in velocity space. 
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1 INTRODUCTION 

' Hierarchical theories of structure formation in the Universe 
propose that galaxies are the result of mergers and accretion 
, of smaller building blocks (White & Rees 1978). Detailed 
' studies of the properties of a galaxy built in this way have 
shown that such events leave fossil signatures in the present 
. day components, which for a galaxy like our own would be 
' clearly detectable with future astrometric missions (Helmi & 
White 1999). In particular the stellar halo would be the nat- 
ural place to look for such substructures, since a spheroidal 
component is formed by the trails of stars left by disrupted 
satellite galaxies. Moreover, recent observations have shown 
that indeed considerable structure is still present in Milky 
Way's halo, indicating that accretion events have had some 
role in its formation history (e.g. Ibata, Gilmore & Irwin 
1994; Majewski, Munn & Hawley 1996; Helmi et al. 1999; 
Ivezic et al. 2000). 

In the next ten years, several satellite missions will be 
devoted to measure with very high accuracy the motions of 
thousands to many millions of stars in our Galaxy. NASA's 
Space Interferometry Mission (SIM) is a targeted mission 
which will obtain parallaxes and proper motions for about 
10000 stars. With somewhat different goals, and more sim- 
ilar to the HIPPARCOS satellite, the Full-sky Astromet- 
ric Mapping Explorer (FAME, Horner et al. 1999) promises 
to measure positions and parallaxes for stars brighter than 
1/ ~ 9 to better than 50 /las and proper motions to 50 /xas 



yr"'^. At V ~ 15 these accuracies will be degraded by an 
order of magnitude. The resulting astrometric database will 
have 4 x 10*^ stars, and may be combined with the radial 
velocities from the Sloan Digital Sky Survey or from other 
ground based catalogues to obtain full phase-space informa- 
tion. Less ambitious but still an improvement over HIPPAR- 
COS is the German DIVA mission (Roser 1998). If launched 
it will observe of the order of 3.5 x 10^ stars, at four times 
the precision of HIPPARCOS (a^ = 0.25 mas and = 0.4 
mas yr~^ ai V = 10), thereby completing the knowledge 
of nearby stars. Like FAME, DIVA will not measure radial 
velocities. On the other hand, the proposed ESA astromet- 
ric satellite GAIA (Gilmore et al. 1998) will provide very 
precise astrometry (<10 pias in parallax and <10 pias yr~^ 
in proper motion at V ~ 15, increasing to 0.2 mas yr~^ at 
V ~ 20) and multicolour photometry, for all 1.3 billion ob- 
jects to V ~ 20, and radial velocities with accuracies of a 
few kms"'^ for most stars brighter than V ~ 17, so that full 
and homogeneous six-dimensional phase-space information 
will be available. These satellite missions will thereby pro- 
vide a very large and statistically reliable sample of stars, 
from which the fundamental questions concerning the origin 
and evolution of the Galaxy may finally be answered. 

In this paper we shall focus on what GAIA will tell 
us about the history and formation of the stellar halo of 
the Milky Way. We will also discuss the impact of DIVA 
and FAME, and leave aside SIM as this mission will not 
provide a survey but a hand-picked catalogue of stars. Even 
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though we focus on the stellar halo, the method that we shall 
propose for finding substructures in phase-space may also be 
extended to find, for example, disk moving groups (e.g. de 
Zeeuw et al. 1999; Chereul, Creze & Bienayme 1999). 

There are several methods for detecting moving groups. 
The Great Circle Cell Counts method (G3C) proposed by 
Johnston, Hernquist & Bolte (1996) uses the position on the 
sky, and employs the fact that satellite galaxies in orbits that 
probe only the outer (spherical) halo conserve the orienta- 
tion of their plane of motion, thereby leaving their debris 
along great circles on the sky, if observed from the Galactic 
centre. The methods used in the Solar neighbourhood for 
detection of disk moving groups and open clusters use also 
proper motions (and sometimes parallax), and assume that 
all the stars belonging to the same system have the same 
velocity vector (e.g. Hoogerwerf & Aguilar 1998; de Bruijne 
1998). Lynden-Bell & Lynden-Bell's method (1995) needs 
the position on the sky and the radial velocity, and has been 
used, for example, to link globular clusters which lie in the 
same plane to some of the (disrupted) dwarf companions of 
our Galaxy (see also Lynden-Bell, 1999). The applicability 
of the above mentioned methods is questionable in the in- 
ner parts of the halo. In this regime, the Galactic potential 
is significantly flattened so that the debris does not remain 
on a fixed plane, the situation where G3C works. As noted 
by Helmi & White (1999) no spatial correlations should be 
expected for satellites disrupted several Gyr ago. On the 
other hand, even though the velocity dispersions in a stellar 
stream do decrease with time, and therefore, very strong 
correlations are to be expected, in the inner halo strong 
phase-mixing takes place. For example in the Solar neigh- 
bourhood several hundred (mainly) cold streams originating 
in disrupted satellites may be present, but it may in prac- 
tice be difficult to resolve each one of such moving groups 
completely. Clearly, before exploring the full capabilities of 
the next generation of astrometric satellite missions, we first 
need to identify where the clustering that is characteristic 
of a satellite manifests itself in the debris that we observe 
after many galactic orbits. As shown in Helmi et al. (1999) a 
method based on the lumpiness in integrals of motion space 
seems to be a promising tool for unveiling the merger history 
of our Galaxy. 



2 BUILDING UP A STELLAR HALO 

Our main goal is to test whether with the next generation 
of astrometric satellites, we would be able to find the signa- 
tures left by merger events in the Galactic stellar halo. We 
will assume that the whole stellar halo is the result of the su- 
perposition of several disrupted satellite galaxies which fell 
onto the Milky Way about 10 Gyr ago. We shall here dis- 
cuss the initial conditions and the numerical methods used 
to generate this version of the stellar halo. 



2.1 Initial conditions for the satellites 

2.1.1 Orbital properties 

The stellar halo has a density profile (Kinman, Suntzeff & 
Kraft 1994) 



P*{r) =Po[ — 
. ro 



(1) 



a total luminosity of about 10 L© , and a half light radius 
which probably lies around 3 kpc from the Galactic centre. 
For r = ro — S kpc (the distance to the Galactic Centre from 
the Sun), po corresponds to the local stellar halo density, for 
which we take po = 1.5 x 10'' kpc"^ (Fuchs & Jahreifi 
1998). 

The initial orbital conditions of our satellites should be 
drawn from the Galactic halo distribution function (DF), 
which we assume to be a function of energy E and angular 
momentum L: f{E,L). For simplicity here we shall assume 
that the stellar halo is a power-law tracer population embed- 
ded in a singular isothermal sphere, representing the dark- 
matter halo of the Milky Way. Following van den Bosch et 
al. (1998), we assume that 



f{E,L)=g{E)h4n), 



where r) = L/Lc{E), 



and Lc(E) is the angular momentum of a circular or- 
bit with energy E: Lc{E) = rc{E)Vc, where rc{E) = 
e'^^"^ exp [iJ/l/i^] . The function ha {rj) is known as the circu- 
larity, and determines the degree of anisotropy of the DF. We 
choose a simple parametrization of hdrf) (Gerhard 1991): 
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so that for a = 0, the DF is isotropic, for a < it is radially 
anisotropic and for a > is is tangentially anisotropic. We 
shall take a = —0.5, since the halo appears to be radially 
anisotropic. 

For a singular isothermal sphere 

The corresponding DF is 
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Here i^max = v&x/us"", with u = (E ~ 4>)/Vc . Since the 
density profile may be derived from the initial distribution 
function as 

r 
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the joint probability distribution of E and r; at a given radius 
rr, is 



(6) 
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(van der Marel, Sigurdsson & Hernquist 1997). Using Eq.(^, 
we find that the normalised cumulative probability distribu- 
tion of E is 



P(< E) = \- exp 



E 



(7) 



We may derive the initial positions of the satellites by 
assuming the profile given in Eq.(^, and using, from Eq.(^, 
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the energy as -E = ~ 14^111(1 — TZ), with TZ a uni- 

form random variable. With the energy (or TZ) we can com- 
pute ryniax = \/2i^(l-7^)^v'-Ml-7^), for 7 = 1/2.5 to 
mimic the stellar halo power law. Using rjuiax and the proba- 
bility distribution for rj we may derive the non-circularity of 
the orbits, and in this way fully determine the phase-space 
initial position of a satellite. 



2.1.2 Internal properties of the satellites 

For the satellites we assume they initially have King pro- 
files, as do most of the satellites in the Local Group. Their 
present day total luminosity is fixed to be that of the stel- 
lar halo. The luminosity of each satellite is drawn from a 
Gaussian distribution with mean 2.5 x 10^ Lq and dispersion 
10'^ L0, thus reproducing the characteristic luminosities of 
Local Group dwarf spheroidals. The initial number of satel- 
lites is thus set to be 33. We assume that the satellites follow 
the scaling relations (Burstein et al. 1997) 

logL = 5.35 + 1.80 log 
log 7? = -0.82 + 0.51 log 

These relations allow us to derive from the luminosity L, 
the core radius R, and the central velocity dispersion a^. 
Our King models have a concentration parameter c = 
logrt/rK ~ 0.72, where n and are the tidal and King 
radii respectively. The initial mass of the satellite is now also 
fully determined. 



2.2 Galactic potentials 

We will consider two different Galactic potentials. In both 
cases, our Galaxy has three components: a dark halo, a disk 
and a bulge, but we take different functional forms for the 
potential. In Model I, we take a dark logarithmic halo 

"l>haio = i^haioln(r^ + (f), (8) 



a Miyamoto-Nagai disk 



"l>dis 



i?2 + (a + Vz^ + b'^Y 



and a spherical Hernquist bulge 



bulge 



r + c 



(9) 



(10) 



where d—12 kpc and fhaio = 131.5 kms ^; Maisk = 10^^ Mq, 
a = 6.5 kpc and b = 0.26 kpc; Mbuigc = 3.4 x 10^° M© and 
c = 0.7 kpc. This choice of parameters gives a flat rotation 
curve with an asymptotic circular velocity of 186 kms~^. 

In Model II, we represent the disk density profile with 
a double exponential (Quinn & Goodman 1986) 



Pd{R,z) = 



Mr 



.g--R/-BDg-/3|z| 



where Rd =3.5 kpc is the disk scale length, Zo is its scale 
height, P = 1/zo and Md = 5.5 x 10^" M© the total disk 
mass. The associated potential is 

GMd 



^d{R, 
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Figure 1. The circular velocity profile as a function of distance 
from the Galactic centre. The solid curve represents the potential 
used in the simulations. The dashed curve corresponds to the 
alternative potential. 
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For the halo we choose (Hernquist 1993) 
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where 

Mh = 1.5 X 10^^ Mq, = (^l-0Fqe«'(l-Erf[g]; 

with q = 'yq/rc and rc = 200 kpc is the cutoff radius. The 
corresponding potential is 



$h(r) = 
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and Ei(2;) is the exponential integral (e. g. G radshteyn & 
Ryzhik 1965). For the bulge we use Eq. (Iicj) but we take 
Mb = 1.1 X 10^° and c = 0.525 kpc (following Velazquez 
& White 1995). Figure 1 shows the circular velocity curves 
produced by each of the two potentials. 

2.3 Numerical methods 

In our numerical simulations, we use the potential of Model I 
for our Galaxy. We represent the satellite galaxy by a col- 
lection of 10^ particles and model their self-gravity by a 
mult ipole expan s ion of the internal pot ential to fourth or- 
der (IWhite 1983] [Zaritsky fc White 1988| ). This type of code 



has the advantage that a large number of particles can be 
followed in a relatively small amount of computer time. In 
this quadrupole expansion, higher than monopole terms are 
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Figure 2. Number counts profile N{r) = r^p{r) for the simu- 
lated stellar halo resulting from the superposition of 33 disrupted 
satellite galaxies, after 12 Gyr of evolution. The straight line rep- 
resents the expected r~^-^ law, arbitrarily shifted. 



softened more strongly. We choose ei ^ 0.2 — 0.25-R for the 
monopole term {R is the core radius of the system) and 
£2 = 2ei for dipole and higher terms and for the centre 
of expansion. The centre of expansion is a particle which, in 
practice, follows the density maximum of the satellite closely 
at all times. 

After letting our satellite relax in isolation, we integrate 
each simulation for ~ 12 Gyr. In Figure]^ we show the final 
particle counts in radial bins N{r) — r'^ p{r) as a function of 
distance from the Galactic centre resulting from the super- 
position of all our experiments. For guidance, we also plot 
the expected r~^'^, arbitrarily shifted. We see that within 
the range of 3 to 30 kpc, our simulations follow relatively 
well the profile. Outside this range we see a sharp drop, due 
to the fact that we are (intentionally) not populating the 
outer halo. Since the properties of the inner stellar halo are 
not so well constrained, we do not worry about the fact that 
we find a shallower slope in the inner few kiloparsecs (this 
is also the result of our initial conditions). More important 
is the fact that our simulations can reproduce very well the 
regime where the astrometric missions promise to give ac- 
curate six dimensional phase-space information. 



2.4 Generating catalogues of halo stars 

To generate an artificial catalogue for the Galaxy we assume 
that each particle in our simulations represents a giant star 
of absolute magnitude Mv = 1. The total number of par- 
ticles in our simulations corresponds well to the expected 
number of giant stars in the Galactic halo (based on the lu- 
minosity function, derived for the age and metallicity char- 
acteristic of halo stars). We prefer to take only giant stars 
at this stage because they are bright enough to be easily 
observable from the Sun. We need to determine a limiting 



Table 1. Estimated precision in parallax (cr^, in /^as) and proper 
motion (tr^ , in fias yr~^) as a function of V magnitude. For FAME 
and DIVA we assume a-,^ = cr^ (based on Horner (1999) and 
Roser (1998), respectively). In the case of GAIA, the estimated 
precisions correspond to a K3 III star with no reddening, and 
increase to 0.2 mas at y ~ 20 (Gilmore et al. 1998). 
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magnitude Viim for our "artificial" catalogue, which we de- 
fine so that all giant stars brighter than Viim have accurate 
full 6-dimensional phase-space information. In the case of 
GAIA, we take Vnm = 15, a limit set by the accuracy in the 
radial velocity. For FAME, Vnm = 12.5 as all stars brighter 
than this magnitude will have relative parallax errors a^/vr 
smaller than (or of the order of) 25%. For DIVA, we take 
Vim = 11 for which a^/ir ~ 0.3. Our GAIA, FAME and 
DIVA catalogues have 386144, 12497 and 1742 "stars" with 
Mv ~ 1 respectively. 

The positions and velocities of each particle are first 
transformed into the observables (a, S, n) and {^a, fJ,s, Vr); 
the expected observational "errors" are then added to the 
parallax, the radial velocity and the proper motion, accord- 
ing to Table ^. For GAIA the precision in the radial veloc- 
ity is taken to be 5 kms~^ for V < 14, and to vary like 
a„ = 1Q{V - 14) + lOkms"^ up to V = 15. Since FAME 
and DIVA will not measure radial velocities on board, for 
these we estimate the error = 15kms~'^, as achievable 
from the ground for such large samples. These "observed" 
quantities are then transformed back to "observed" positions 
and velocities. We repeat this procedure 5 times to obtain 5 
different realizations of the data. 



3 FINDING DISRUPTED GALAXIES 
3.1 Integrals of motion space 

Our satellites disrupt relatively quickly, in only a few peri- 
centric passages. Therefore we may consider each of the 33 
satellites as an ensemble of particles with very similar inte- 
grals of motion (energy, angular momentum). As we show 
in Figure ^, initially satellites are both clumps in configura- 
tion and velocity space, as they are in {E, L, Lz) space. If 
these are conserved quantities, or evolve only slightly, this 
initial clumping should be present even after the system has 
phase-mixed completely. Thus the space of integrals or adi- 
abatic invariants seems to be the natural space to look for 
the substructure produced by an accreted satellite. 

There are a few issues we should address here before 
fully discussing a method based on clumping in the integrals 
of motion space. To compute the energy of the particles (or 
stars that will be observed by GAIA for example) we need 
to assume a Galactic potential. To determine the success 
of such a method we need to understand how our lack of 
knowledge on the precise form of the Galactic potential in- 
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Figure 5. Final distribution of particles in the Lz - E space 
after error convolution for FAME (left panel) and for DIVA 
(right panel), with energies computed using the original poten- 
tial. A comparison to the left panel of Figs. ^ and ^shows that 
the expected errors for these missions tend to erase much of the 
substructure left in the integrals of motion space. 
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Figure 6. Final distribution of particles in the Lz - E space af- 
ter GAIA error convolution for the alternative potential. Com- 
pare to left panel of Fig.M. 
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fluences our results. We shall therefore proceed in two steps. 
In the first step, we take the same potential as that used in 
the simulations, Model I. In the final step we use the alter- 
native potential introduced in Sec. 2.2, our Model II. This 



last step, in which we do not know the exact form of the 
Galactic potential but we make a reasonable guess, is most 
likely to represent the real situation. 

Secondly, even though the total angular momentum is 
not fully conserved for an axisymmetric potential (only Lz 
is), it evolves preserving a certain degree of coherence. The 
advantage of using the integrals of motion space is that the 
number of clumps detected in this way will represent well the 
total number of accretion/merging events, since unlike other 
methods which are only local, it singles out all the stars from 
a given accreted object, independently of how different their 
phases and velocities might be. We choose to make use of 
all three integrals to reduce the chances of overlap amongst 
different lumps, since this probability clearly depends on the 
dimensionality of the space. 

The analogue of Figure ^ for particles "brighter than 
15th magnitude" (roughly within 6 kpc from the Sun) in 
the simulations, after 12 Gyr of evolution and for the origi- 
nal potential, shows that, even though there is some degree 
of evolution, clurnping remains in the integrals of motion 
space. In Figure Ml we plot the integrals of motion space 
for one realization of the GAIA catalogue, i.e. after error 
convolution. A number of substructures are clearly visible, 
many of which can be directly related to the initial distribu- 
tion, even with the GAIA observational uncertainties taken 
into account. This shows that the expected observational er- 
rors for GAIA will not affect the chances of detecting such 
substructures. In the case of FAME the situation is not as 
good, as illustrated in the left panel of Figure ^ where the 
different lumps are less populated (because of the magni- 
tude limit) and considerably more smeared out (because of 
the larger observational errors). For DIVA the clumping has 
disappeared almost completely, as shown in the right panel 
of the same figure. 

Figure ^ corresponds to the same realization of the 
GAIA catalogue as used before, but with the energies cal- 
culated using the case of the potential of Model II. Clearly, 
even though the two considered potentials are different, the 
substructure remains. The uncertainty in the precise form 
of the Galactic potential therefore does not affect the likeli- 
hood of finding disrupted satellites. 

3.2 Method: EOF in integrals of motion space 

We use a Friends-of-Friends (FOF) algorithm to find clumps 
in the integrals of motion space. This method has been used 
frequently to find bound halos in cosmological A'^-body simu- 
lations. The basic idea is that all particle pairs separated by 
less than a fraction i of the mean interparticle distance are 
linked. Disjoint sets of connected particles are then identified 
as halos (Efstathiou et al. 1988) . These halos correspond ap- 
proximately to the regions interior to isodensity contours at 
an overdensity of 2/£^. This FOF procedure allows a rapid 
identification of halos, and moreover, all members of a given 
halo found for a particular value of £ are members of the 
same halo in any list generated for a larger value of ^. In 
the case of cosmological simulations, the linking distance is 
defined so that the mean density of a halo is about 200 times 
the density of the Universe at the time of identification. 



In our case, it is less clear how we should define the 
interparticle distance, or the linking length. Because the en- 
ergy and angular momentum have, by definition, different 
scales, it seems natural to try to reduce everything to the 
same scale, or equivalently, to use instead of spheres an el- 
lipsoidal configuration. Even though the angular momentum 
and its z-component have the same scale, lumps are gener- 
ally elongated in the L-direction with a 2:1 ratio, as can be 
seen from Figure |^. We therefore search for lumps whose 
characteristic size would be defined as: 



AL ~ 2ALz 



AE 



(kms-i)2 



20 



kpckms^i ' 



where now AL, would be related to the linking length. This 
implies that we re-scale the variables according to 



E 



E/20, L L/2, L, 



The factor 20 in the energy scaling may be derived (heuristi- 
cally) from the fact that the typical energy range in the Solar 
neighbourhood is 1.6 x 10^(kms^^)^, whereas the range of 
Lz is 8000 kpc kms"^ (from -4000 to 4000 kpc kms"^). 

We will apply the FOF algorithm for two different link- 
ing lengths, to allow for different characteristic sizes of the 
halos and resolutions in the algorithm. Note that there are 
particular regions in this space which are occupied by more 
than one satellite, even in this 3-dimensional space (this is 
even worse if only the Lz — L plane is used) , so that not each 
of the lumps found may correspond to only one satellite, but 
may have contributions of a few. 

3.3 Results 

We apply the FOF algorithm to our GAIA catalogue, in- 
cluding error convolution for all particles brighter than 15th 
magnitude, and using the original potential. We take two 
different values for the FOF linking length: ^ = 16 and 
I — 30, where AL^ — 51. We consider groups with at least 
500 "stars" . We combine the two group catalogues to obtain 
a new group catalogue which contains all lumps detected. If 
some particles are found to belong to two different clumps 
(one from each catalogue) we keep the lump which has the 
smallest size. We now iterate one more time on the cata- 
logue defined by the particles that do not belong to any 
of the lumps found by our FOF, again for the two values 
of £ and with a minimum of 250 "stars" . Some of the newly 
found lumps can be related to those previously detected, and 
some others are found to resolve some of the largest lumps 
in our initial group catalogue. In the left panel of Figure ^ 
we show the distribution of energy E and L, for our final 
group catalogue. 

We find 17 different groups with this method. Not all 
the groups may be associated exclusively with one of our 
original satellites. As can be seen from Figj^, there is quite 
a bit of superposition in this three-dimensional space, and 
so not all the original satellites can be recovered, or equiv- 
alently, not all lumps can be resolved with just two itera- 
tions. If we analyse how the particles in the different lumps 
can be related to particles in the initial satellites we find 
that, out of the 17 groups discovered, 14 can be associated 
almost uniquely to one satellite^. This means that our sim- 



We say that a group is almost uniquely associated to one satel- 
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Figure 7. Lumps detected with our FOF algorithm. In panel (a) we show the final group catalogue for the original potential used in 
the simulations after convolving with the observational errors expected for GAIA. Panel (b) corresponds to our alternative potential and 
also to the GAIA catalogue. In both cases the recovery rate is about 50%. Panel (c) shows the lumps recovered by our FOF apnlicd to 
the FAME catalogue generated as described in the text and for energies computed with the original potential. Compare to Fig. Win the 
case of GAIA and to the left panel in Fig.0 for FAME. 
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Figure 8. The velocity space distribution for particles in a cubic volume of 2 kpc on a side centered on the Sun, and for one realization 
of the GAIA catalogue. In the upper panels different colours indicate particles associated with different satellites (using the same colour 
coding as in Figure H). In the lower panels, the colours are used to show particles associated to the lumps recovered by our FOF algorithm 
applied to the GAIA catalogue in the case of the alternative Galactic potential. (Here the colour coding corresponds to that used in 
panel (b) of Figure H.) 
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pie method is capable of finding more than 40% of all the 
satellites that were accreted by our "Galaxy" . 

Similarly, we apply the FOF algorithm to the same 
GAIA catalogue but now compute the energy E of the par- 
ticles with the alternative potential. In this case we again 
find 17 different lumps (after two iterations, and combining 
the results of the two different values of the linking length). 
Of these 17 groups, 14 can be uniquely associated to one 
satellite. This is shown in the central panel of Figure ^ 
Our method is thus quite successful in identifying disrupted 
satellite galaxies in integrals of motion space, even when we 
only have a guess for the Galactic potential and when the 
observational uncertainties are taken into account. Further 
uncertainties such as the distance from the Sun to the Galac- 
tic centre and the velocity of the Local Standard of Rest are 
also negligible. 

When we apply the same method for the original poten- 
tial on the FAME catalogue we are able to find 6 different 
groups. For 5 of these a unique correspondence with an ac- 
creted satellite exists, as shown in the rightmost panel of 
Figure ^. In the case of DIVA we find only 1 group (with 
at least 20 particles), which can be easily identified visu- 
ally from the right panel in Figure In the cases of DIVA 
and FAME we used slightly larger linking lengths to take 
into account the smearing out of the lumps caused by the 
larger observational uncertainties. Because the samples are 
also smaller (because of the limiting magnitude), we con- 
sider groups with at least 50 particles in the case of FAME, 
and 20 particles for DIVA. 



4 CLUMPINESS IN THE KINEMATICS OF 
HALO STARS 

In Figure ^ we show the velocities of all the particles con- 
tained in a volume of 2 kpc on a side for one realization 
of the GAIA catalogue. There is considerable substructure, 
which is visible thanks to the great precision that GAIA 
will achieve. From the upper panels it is clear that, as dis- 
cussed in the introduction, distinguishing the satellites that 
gave rise to each one of the different moving groups is a 
non-trivial task in this space. In the lower panels we have 
coloured the different contributions from the 14 groups de- 
tected by our FOF algorithm in the case of the alternative 
potential. A comparison between upper and lower panels 
also shows how successful our method is. 

The kinematically cold streams visible in Figure ^ re- 
main as coherent structures for longer than a Hubble time. 
This is true even when mergers, rather than simple satellite 
accretion, are dominant (Helmi, White & Springel 2000). 
The dumpiness in the kinematics of halo stars should thus 
be a distinct feature of the hierarchical formation of our 
Galaxy. It is therefore also interesting to determine the de- 
gree of the dumpiness and whether it will be measurable 
with future astrometric missions. We will determine this 
dumpiness using the two-point correlation function ^ in ve- 
locity space for a sphere of 1 kpc radius around the Sun. We 
estimate ^ from 



lite if more than 70% of the particles in the group belong to only 
that satellite. 



^_ {DD){RR) 



(14) 



(e.g. Hamilton 1993) where {DD} is the normalised number 
of pairs of particles with velocities in a given velocity range 
(or bin), i.e. 



(DD) 



pairs of particles i,j with li < | V; — Vj | < f -|- A 
" NoiNo - 1) 



(15) 



and where Nd is the number of particles in the sphere. (RR) 
is defined analogously but for Nr random points. The ran- 
dom variates are drawn from a trivariate Gaussian distri- 
bution determined from the "data" in the principal axes 
velocity frame. Here we take Nr = IONd ■ Finally (DR) are 
the normalised counts for "data" -random pairs. We estimate 
the uncertainty from 



Ae = (l + C) 



Nd{Nd-1){DD)' 



(16) 



If the sample contains kinematically cold streams, we 
should find an excess of pairs in the bins corresponding to 
small velocity differences, i.e. the correlation function should 
be significantly different from zero (which corresponds to the 
absence of correlations). We proceed by measuring ^ for our 
GAIA, FAME and DIV A ca talogues including error convolu- 
tion as described in Sec. 



2.4 



We also vary the position of the 
1 kpc sphere around the "Sun", keeping the same distance 
from the "Galactic centre". That is, we place the Sun at 
{x,y) = {(8,0),(0, -8),(-8,0),(0,8)kpc} and z = 0. This 
allows us to account for the natural variations one may have 
from volume to volume. We then make five realizations for 
each "Sun" position for each catalogue. In Figure^ we show 
the correlation function obtained by averaging over all the 
realizations, for each 1 kpc sphere and for each catalogue. 
The average ^ for each volume is the weighted mean, where 
the weights are given by 1/A|, and the error bars indicate 
the (weighted) dispersion around the (weighted) mean. We 
find an excess of pairs of stars with similar motions, the sig- 
nature indicating the presence of cold streams as expected 
for a stellar halo built by disrupted satellites. Note that it 
will even be possible to determine that the halo is not a 
smooth distribution in the Solar neighbourhood even with 
velocity errors of the order of 20 kms~^ such as those ex- 
pected for DIVA for a star with Mv = 1 at 500 pc from the 
Sun. 



5 DISCUSSION 

We simulated the entire stellar halo of the Galaxy start- 
ing from disrupted satellite galaxies, and "observed" it with 
the next generation of astrometric satellites (DIVA, FAME 
and GAIA). We analysed the observations with the aim of 
recovering the different accretion events our "Galaxy" expe- 
rienced over its lifetime. We used a FOF algorithm to find 
clumps in the integrals of motion space, which we expected 
would correspond to the disrupted satellites. Our integrals 
of motion space is defined by energy E, total angular mo- 
mentum L and its z-component Lz, even though strictly 
speaking these are not fully conserved quantities (because 
of interaction of the stars while still bound to the satellite, 
and because of the axisymmetry of our Galaxy). We have 
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Figure 9. The two-point correlation function for "giant stars" inside spheres of 1 kpc radius around the Sun (defined as 8 kpc from the 
Galactic centre on the Galactic disk) computed as the weighted average over five realizations of the DIVA, FAME and GAIA catalogues. 
The different symbols correspond to f measured inside spheres at different locations of the "Sun" on the Solar circle. 



shown that the initial clumping in this space is maintained 
to a great extent even after 12 Gyr of evolution. 

After using our FOF algorithm we find that we can only 
recover a couple of accreted satellites (in our analysis just 
one) for the DIVA catalogue, whereas for FAME we recover 
about 15% of all satellites. In both cases we assume that the 
astrometry is complemented by ground based radial velocity 
measurements. The situation is significantly different in the 
case of the GAIA catalogue, for which we recover almost half 
of all disrupted satellites with this simple algorithm. The im- 
provement generally lies in the larger volume for which full 
6-D information is available, in particular when comparing 
FAME and GAIA. The use of 6-D information appears to be 
essential to recover all the events, as there is a large fraction 
of phase-space where these are superposed. This is particu- 
larly clear from Figure ^ (rightmost panel), where angular 
momentum alone cannot be used to distinguish the different 
satellites. Whereas by eye inspection in the {E, L, L^) space 
we may recover five or six events, for the space (L, Lz) this is 
reduced to one or two events. Our results are unlikely to be 
strongly dependent on the particular choice of the luminos- 
ity distribution of the disrupted satellites. This is because 
a large number of small satellites occupy basically the same 
phase-space volume as a small number of large ones. 

The evolution of the Galactic potential may be the most 
crucial simplification in our analysis. In hierarchical cos- 
mologies the number of objects that form a galaxy like our 
own is in the range of 5 — 20, with comparable masses. The 
process of formation is likely to be very violent and the po- 
tential is surely not static, quite probably not axisymmet- 
ric, and therefore the initial clumping of the system may 
not be reflected in clumping in our defined integrals of mo- 
tion space. However, if this happened during the first few 



Gyrs, any object infalling later, ought to have perceived a 
fairly static (or adiabatically changing) Galaxy, and then 
our method would still be useful. Indeed, some preliminary 
analysis of the formation of a halo in a ACDM cosmology 
indicates that particles from different satellites may be re- 
covered as lumps in this space (Helmi et al. 2000), though 
the structure is less evident than in the plots shown here, 
where even by simple eye inspection one may recover about 
1/5 of all satellites. 

What will anyway remain as signatures of the merger 
history of our Galaxy will be the kinematically cold streams 
originating in disrupted halos. An interesting observational 
test is the comparison of the kinematics of a smooth, pos- 
sibly Gaussian, distribution (which may be expected in the 
case of a monolithic collapse) to the kinematics observed in 
the stellar halo built by disrupted satellites. Our analysis of 
the correlation function in velocity space indicates the pres- 
ence of a larger number of streams with very small velocity 
dispersions in a sphere of 1 kpc radius around the Sun. This 
test will be feasible even for DIVA. The key to the success of 
this test lies in the complete and large sample of stars with 
3-D velocities which will be available. 

In this paper we have focused on determining the 
merger history of the Milky Way, rather than the precise 
form of the Galactic potential or to what extent it may have 
varied. However, these are key questions that will be solved 
very likely by SIM and GAIA (e.g. Johnston et al. 1999). 
We may add here that after finding the different satellites 
we will be able to determine the conditions and character- 
istics of objects that fell onto the Milky Way more than 10 
Gyr ago. 
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